In [47]:
%reload_ext autotime
import rasterio
import requests
from rasterio.features import shapes, sieve, rasterize
from rasterio.plot import show
import geopandas as gpd
import os
import numpy as np
from tqdm.auto import tqdm
from shapely.geometry import box, shape
from rasterio.mask import mask
from matplotlib import pyplot as plt
import xdem
import pandas as pd
tqdm.pandas()
pd.set_option('display.max_columns', None)
✔️ 2.4 ms (2024-12-02T08:53:27/2024-12-02T08:53:27)
In [109]:
manifest = pd.json_normalize(requests.get("https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/collection.json").json()["links"])
manifest = manifest[manifest["rel"] == "item"]
manifest["tilename"] = manifest.href.str.replace(".json", "").str.strip("./")
manifest
✔️ 233 ms (2024-12-02T10:00:26/2024-12-02T10:00:26)
Out[109]:
| rel | href | type | file:checksum | tilename | |
|---|---|---|---|---|---|
| 2 | item | ./BD43_10000_0204.json | application/json | 12208b52bfdff84d924430e9503bc70bcc7555c1861e5b... | BD43_10000_0204 |
| 3 | item | ./BD43_10000_0205.json | application/json | 12204ec4b1257e965c8da142c4696b1f8cbbc6c78d0cf0... | BD43_10000_0205 |
| 4 | item | ./BD43_10000_0304.json | application/json | 1220acfd4fd86a3866f4facf92bb75cd5224451e774ff4... | BD43_10000_0304 |
| 5 | item | ./BD43_10000_0305.json | application/json | 122051507adf6a042f09c138f1842911fec7bd43c270d8... | BD43_10000_0305 |
| 6 | item | ./BD43_10000_0404.json | application/json | 1220c63b7011af35652df009f8a6b69f4e3f8ea2db33d6... | BD43_10000_0404 |
| ... | ... | ... | ... | ... | ... |
| 302 | item | ./BH43_10000_0102.json | application/json | 12200f2d3649a677e745f258162ca6ebc13904d6ecde95... | BH43_10000_0102 |
| 303 | item | ./BH43_10000_0201.json | application/json | 1220d28f9ff3b2f6ad76139b4ccfa5bd86bc996b417e36... | BH43_10000_0201 |
| 304 | item | ./BH43_10000_0202.json | application/json | 1220d78d97c89b75b5d35667f77acb8e5a2078d52b977b... | BH43_10000_0202 |
| 305 | item | ./BH43_10000_0301.json | application/json | 12200b058bab70e426eccd9ae990a5888e1754b27b97c2... | BH43_10000_0301 |
| 306 | item | ./BH43_10000_0302.json | application/json | 12203c93b413b89941df279575d9d352046cc8b6977deb... | BH43_10000_0302 |
305 rows × 5 columns
In [2]:
if not os.path.isfile("old.tif"):
with open("old.tif", "wb") as f:
response = requests.get("https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2018-2020/dem_1m/2193/BD44_10000_0503.tiff")
f.write(response.content)
with open("new.tif", "wb") as f:
response = requests.get("https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/BD44_10000_0503.tiff")
f.write(response.content)
old = rasterio.open("old.tif")
new = rasterio.open("new.tif")
✔️ 23.8 ms (2024-11-28T14:29:32/2024-11-28T14:29:32)
In [3]:
new.crs, old.crs
✔️ 5.51 ms (2024-11-28T14:29:32/2024-11-28T14:29:32)
Out[3]:
(CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'),
CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'))
In [4]:
show(old, title="2018-2020 DEM")
✔️ 2.07 s (2024-11-28T14:29:32/2024-11-28T14:29:34)
Out[4]:
<Axes: title={'center': '2018-2020 DEM'}>
In [5]:
show(new, title="2023 DEM")
✔️ 1.72 s (2024-11-28T14:29:34/2024-11-28T14:29:36)
Out[5]:
<Axes: title={'center': '2023 DEM'}>
In [6]:
diff = new.read(1) - old.read(1)
print(diff.shape)
plt.imshow(diff, cmap="coolwarm_r", vmin=-1, vmax=1)
plt.colorbar()
✔️ 1.54 s (2024-11-28T14:29:36/2024-11-28T14:29:37)
(7200, 4800)
Out[6]:
<matplotlib.colorbar.Colorbar at 0x7f4453e3cd00>
In [7]:
_ = plt.hist(diff.flatten(), bins=100, range=(-2, 2))
plt.title("Histogram of elevation differences")
✔️ 516 ms (2024-11-28T14:29:37/2024-11-28T14:29:38)
Out[7]:
Text(0.5, 1.0, 'Histogram of elevation differences')
In [8]:
result = diff.round().clip(min=-1, max=1).astype(np.int16)
result = sieve(result, 4000)
result = np.where(result != 0, result, np.nan)
plt.imshow(result, cmap="coolwarm_r")
plt.colorbar()
✔️ 3.41 s (2024-11-28T14:29:38/2024-11-28T14:29:41)
Out[8]:
<matplotlib.colorbar.Colorbar at 0x7f4453f6f820>
In [9]:
areas = gpd.GeoDataFrame(geometry=[shape(s) for s, v in shapes(result, transform=new.transform) if not np.isnan(v)])
areas["area"] = areas.area
areas.sort_values("area", ascending=False, inplace=True)
areas.plot("area", legend=True)
✔️ 1.28 s (2024-11-28T14:29:42/2024-11-28T14:29:43)
Out[9]:
<Axes: >
In [53]:
largest_area = areas.iloc[0]#.head(1)
largest_area
✔️ 10.6 ms (2024-12-02T09:01:57/2024-12-02T09:01:57)
Out[53]:
geometry POLYGON ((2054419 5803910, 2054419 5803909, 20... area 530308.0 Name: 176, dtype: object
In [55]:
masked_old, transform = mask(old, [largest_area.geometry], nodata=np.nan)
nonan = np.argwhere(~np.isnan(masked_old[0]))
top_left = nonan.min(axis=0)
bottom_right = nonan.max(axis=0)
print(top_left, bottom_right)
masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_old)
✔️ 1.33 s (2024-12-02T09:02:28/2024-12-02T09:02:30)
[5290 0] [7044 3057]
Out[55]:
<Axes: >
In [56]:
masked_new, transform = mask(new, [largest_area.geometry], nodata=np.nan)
masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_new)
✔️ 1.1 s (2024-12-02T09:02:38/2024-12-02T09:02:40)
Out[56]:
<Axes: >
In [26]:
features = {
"old_min": np.nanmin(masked_old),
"old_max": np.nanmax(masked_old),
"old_mean": np.nanmean(masked_old),
"old_median": np.nanmedian(masked_old),
"old_std": np.nanstd(masked_old),
"new_min": np.nanmin(masked_new),
"new_max": np.nanmax(masked_new),
"new_mean": np.nanmean(masked_new),
"new_median": np.nanmedian(masked_new),
"new_std": np.nanstd(masked_new),
}
features
✔️ 149 ms (2024-11-28T14:35:27/2024-11-28T14:35:27)
Out[26]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585}
In [22]:
raster = rasterize(largest_area.geometry, out_shape=result.shape, transform=new.transform)
masked_diff = np.where(raster, diff, np.nan)
masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_diff)
✔️ 617 ms (2024-11-28T14:33:52/2024-11-28T14:33:52)
Out[22]:
<Axes: >
In [27]:
features.update({
"diff_min": np.nanmin(masked_diff),
"diff_max": np.nanmax(masked_diff),
"diff_mean": np.nanmean(masked_diff),
"diff_median": np.nanmedian(masked_diff),
"diff_std": np.nanstd(masked_diff),
})
features
✔️ 76 ms (2024-11-28T14:35:30/2024-11-28T14:35:30)
Out[27]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585,
'diff_min': -6.8740234,
'diff_max': 43.521973,
'diff_mean': 8.370939,
'diff_median': 4.911972,
'diff_std': 9.06017}
In [20]:
attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
attributes = xdem.terrain.get_terrain_attribute(
masked_diff,
resolution=old.res,
attribute=attribute_names
)
plt.figure(figsize=(8, 6.5))
for i in range(8):
plt.subplot(4, 2, i + 1)
plt.imshow(attributes[i].squeeze(), cmap="viridis")
cbar = plt.colorbar()
cbar.set_label(attribute_names[i])
plt.xticks([])
plt.yticks([])
plt.tight_layout()
plt.show()
✔️ 20.4 s (2024-11-28T14:29:55/2024-11-28T14:30:15)
In [29]:
for i, name in enumerate(attribute_names):
features.update({
f"{name}_min": np.nanmin(attributes[i]),
f"{name}_max": np.nanmax(attributes[i]),
f"{name}_mean": np.nanmean(attributes[i]),
f"{name}_median": np.nanmedian(attributes[i]),
f"{name}_std": np.nanstd(attributes[i]),
})
features
✔️ 841 ms (2024-11-28T14:36:17/2024-11-28T14:36:18)
Out[29]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585,
'diff_min': -6.8740234,
'diff_max': 43.521973,
'diff_mean': 8.370939,
'diff_median': 4.911972,
'diff_std': 9.06017,
'roughness_min': 0.024993896484375,
'roughness_max': 16.99798583984375,
'roughness_mean': 0.8664416971654487,
'roughness_median': 0.501007080078125,
'roughness_std': 0.9008248692669407,
'slope_min': 0.00030909907704678905,
'slope_max': 83.09449329265449,
'slope_mean': 15.825287450127606,
'slope_median': 9.967977228482425,
'slope_std': 14.601294984565905,
'aspect_min': 0.0,
'aspect_max': 359.9997588896344,
'aspect_mean': 172.74772778897523,
'aspect_median': 183.00927436545985,
'aspect_std': 104.09906838438863,
'curvature_min': -1643.609619140625,
'curvature_max': 2639.605712890625,
'curvature_mean': 1.7315227064581844,
'curvature_median': 0.8087158203125,
'curvature_std': 43.72569863077058,
'terrain_ruggedness_index_min': 0.024669455802839076,
'terrain_ruggedness_index_max': 27.121265279417575,
'terrain_ruggedness_index_mean': 0.8939862986187255,
'terrain_ruggedness_index_median': 0.5406418612709043,
'terrain_ruggedness_index_std': 0.8938326265050937,
'rugosity_min': 1.000078699472661,
'rugosity_max': 9.8867222334156,
'rugosity_mean': 1.1083063933475281,
'rugosity_median': 1.0261556432873709,
'rugosity_std': 0.19958467669264315,
'profile_curvature_min': -1767.2495203822261,
'profile_curvature_max': 1433.1465250703304,
'profile_curvature_mean': -1.2548662332856713,
'profile_curvature_median': -0.5842965641336939,
'profile_curvature_std': 28.49170269374419,
'planform_curvature_min': -624.432616658446,
'planform_curvature_max': 978.0291066330788,
'planform_curvature_mean': 0.4766817029082818,
'planform_curvature_median': 0.18651537440086147,
'planform_curvature_std': 21.37279648106214}
In [82]:
REC = gpd.read_file("nzRec2_v5.gdb", engine='pyogrio', use_arrow=True)
REC
✔️ 8.05 s (2024-12-02T09:39:23/2024-12-02T09:39:31)
/usr/local/lib/python3.10/dist-packages/pyogrio/raw.py:287: UserWarning: More than one layer found in 'nzRec2_v5.gdb': 'riverlines' (default), 'Hydro_Net_Junctions', 'rec2ws', 'rec2_rain_runoff_V5', 'rec2_stew_baserock_V5', 'rec2_stew_toprock_V5', 'rec2_utility_variables_V5', 'rec2_lcdb2_V5', 'rec2_lcdb1_V5', 'rec2_lcdb3_V5', 'rec2_ni_baserock_V5', 'rec2_ni_toprock_V5', 'rec2_si_baserock_V5', 'rec2_si_toprock_V5', 'rec2_elevband_rain_V5', 'LAYERKEYTABLE', 'APUNIQUEID', 'riverlines_FS', 'N_1_Desc', 'N_1_EDesc', 'N_1_EStatus', 'N_1_ETopo', 'N_1_FloDir', 'N_1_JDesc', 'N_1_JStatus', 'N_1_JTopo', 'N_1_JTopo2', 'N_1_Props'. Specify layer parameter to avoid this warning. with open_arrow(
Out[82]:
| OBJECTID | HydroID | NextDownID | CATAREA | CUM_AREA | nzsegment | Enabled | LENGTHDOWN | Headwater | Hydseq | StreamOrde | euclid_dis | upElev | downElev | upcoordX | downcoordX | downcoordY | upcoordY | sinuosity | nzreach_re | headw_dist | segslpmean | LID | reachtype | FROM_NODE | TO_NODE | Shape_Leng | Shape_Length | FLOWDIR | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 9 | 218553.786394 | 2.185538e+05 | 1000005 | 1 | 1801.260253 | 1 | 1 | 1 | 199.150697 | 119.994514 | 90.628487 | 1.601529e+06 | 1.601574e+06 | 6.192386e+06 | 6.192580e+06 | 1.072393 | 1000004 | 0 | 8.388180 | 0 | 0 | 1 | 2 | 213.567866 | 213.567866 | 1 | MULTILINESTRING ((1601528.857 6192580.400, 160... |
| 1 | 2 | 2 | 9 | 455995.848576 | 4.559958e+05 | 1000003 | 1 | 1801.260253 | 1 | 2 | 1 | 537.313689 | 163.071075 | 90.628487 | 1.601783e+06 | 1.601574e+06 | 6.192386e+06 | 6.192881e+06 | 1.082775 | 1000005 | 0 | 7.678527 | 0 | 0 | 3 | 2 | 581.789877 | 581.789877 | 1 | MULTILINESTRING ((1601782.856 6192881.076, 160... |
| 2 | 3 | 3 | 4 | 461385.667778 | 4.613857e+05 | 1000008 | 1 | 335.656220 | 1 | 4 | 1 | 596.782205 | 64.632355 | 20.004650 | 1.599355e+06 | 1.599237e+06 | 6.191779e+06 | 6.192364e+06 | 1.163688 | 1000009 | 0 | 4.276652 | 0 | 0 | 4 | 5 | 694.468045 | 694.468045 | 1 | MULTILINESTRING ((1599355.244 6192363.831, 159... |
| 3 | 4 | 4 | -1 | 135806.633398 | 2.330307e+06 | 1000010 | 1 | 0.000000 | 0 | 9 | 2 | 288.064229 | 20.004650 | 4.172440 | 1.599237e+06 | 1.598982e+06 | 6.191913e+06 | 6.191779e+06 | 1.165213 | 1000010 | 788 | 3.145852 | 0 | 0 | 5 | 6 | 335.656220 | 335.656220 | 1 | MULTILINESTRING ((1599237.075 6191778.666, 159... |
| 4 | 5 | 5 | -1 | 863421.445984 | 8.634214e+05 | 1000006 | 1 | 0.000000 | 1 | 5 | 1 | 1020.122542 | 143.468826 | 4.437851 | 1.602684e+06 | 1.602267e+06 | 6.191623e+06 | 6.192554e+06 | 1.182202 | 1000006 | 0 | 7.760943 | 0 | 0 | 7 | 8 | 1205.990464 | 1205.990464 | 1 | MULTILINESTRING ((1602683.557 6192553.928, 160... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 593512 | 593614 | 593513 | -1 | 51990.267885 | 1.174888e+09 | 4170001 | 1 | NaN | 0 | 124017 | 6 | 0.000000 | 0.000000 | 0.000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 85116 | 604841 | 542.278140 | 542.278140 | 1 | MULTILINESTRING ((1975828.498 5785372.292, 197... |
| 593513 | 593616 | 593514 | -1 | 150666.238187 | 2.061735e+08 | 1040001 | 1 | NaN | 0 | 27152 | 5 | 0.000000 | 12.272364 | 0.000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 26864 | 604842 | 1675.015974 | 1675.015974 | 1 | MULTILINESTRING ((1731848.082 6017486.385, 173... |
| 593514 | 593620 | 593515 | -1 | 53587.929573 | 5.503258e+07 | 7260002 | 1 | NaN | 0 | 239462 | 5 | 0.000000 | 4.917045 | 3.455193 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 244056 | 604843 | 156.868127 | 156.868127 | 1 | MULTILINESTRING ((1789499.342 5528919.920, 178... |
| 593515 | 593621 | 593516 | 593515 | 8667.469549 | 5.432822e+07 | 7260001 | 1 | NaN | 0 | 239461 | 5 | 0.000000 | 6.187159 | 4.917045 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 244073 | 244056 | 116.148688 | 116.148688 | 1 | MULTILINESTRING ((1789607.824 5528878.421, 178... |
| 593516 | 593622 | 593517 | -1 | 492400.109562 | 1.856727e+08 | 7247115 | 1 | 529.317107 | 0 | 249508 | 5 | 964.116694 | 5.069054 | 1.886685 | 1.782739e+06 | 1.782650e+06 | 5.496834e+06 | 5.497794e+06 | 1.086063 | 7047051 | 36796 | 0.272352 | 0 | 0 | 252501 | 604844 | 151.702564 | 151.702564 | 1 | MULTILINESTRING ((1782739.481 5497793.710, 178... |
593517 rows × 30 columns
In [84]:
LCDB = gpd.read_file("lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gpkg", engine='pyogrio', use_arrow=True)
LCDB
✔️ 21.6 s (2024-12-02T09:40:19/2024-12-02T09:40:41)
Out[84]:
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | Wetland_18 | Wetland_12 | Wetland_08 | Wetland_01 | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb2000096676 | MULTIPOLYGON (((1613722.435 5425797.372, 16137... |
| 1 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Regional Council | 2019-12-01 | lcdb1000513359 | MULTIPOLYGON (((1816770.219 5947804.627, 18167... |
| 2 | Mangrove | Mangrove | Mangrove | Mangrove | Mangrove | 70 | 70 | 70 | 70 | 70 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb1000182160 | MULTIPOLYGON (((1715672.186 5958842.706, 17156... |
| 3 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb1000065930 | MULTIPOLYGON (((1705330.918 6088979.740, 17053... |
| 4 | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | 22 | 22 | 22 | 22 | 22 | no | no | no | no | no | no | no | no | no | no | Regional Council | 2019-12-01 | lcdb1000065472 | MULTIPOLYGON (((1761684.636 5789742.527, 17616... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 511099 | Low Producing Grassland | Low Producing Grassland | Low Producing Grassland | Manuka and/or Kanuka | Low Producing Grassland | 41 | 41 | 41 | 52 | 41 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb1000505810 | MULTIPOLYGON (((1785112.194 5684560.595, 17851... |
| 511100 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb2000193885 | MULTIPOLYGON (((1607510.750 5432591.699, 16075... |
| 511101 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb2000219027 | MULTIPOLYGON (((1603592.310 5269947.382, 16036... |
| 511102 | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | 45 | 45 | 45 | 45 | 45 | yes | yes | yes | yes | yes | yes | yes | yes | yes | yes | Regional Council | 2014-06-30 | lcdb1000417326 | MULTIPOLYGON (((1822629.596 5477414.337, 18226... |
| 511103 | High Producing Exotic Grassland | High Producing Exotic Grassland | Exotic Forest | Exotic Forest | Exotic Forest | 40 | 40 | 71 | 71 | 71 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30 | lcdb1000145507 | MULTIPOLYGON (((1704402.741 5624819.729, 17044... |
511104 rows × 24 columns
In [91]:
LCDB[LCDB.intersects(largest_area.geometry)]
✔️ 419 ms (2024-12-02T09:46:00/2024-12-02T09:46:00)
Out[91]:
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | Wetland_18 | Wetland_12 | Wetland_08 | Wetland_01 | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 17220 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | High Producing Exotic Grassland | 71 | 71 | 71 | 71 | 40 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000053082 | MULTIPOLYGON (((2054973.343 5805787.881, 20549... |
| 18662 | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | 16 | 16 | 16 | 16 | 16 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30 | lcdb1000016122 | MULTIPOLYGON (((2048258.705 5805949.487, 20482... |
| 20295 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127897 | MULTIPOLYGON (((2054496.624 5803468.778, 20545... |
| 23069 | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | 54 | 54 | 54 | 54 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb1000119141 | MULTIPOLYGON (((2054140.516 5804625.890, 20541... |
| 349500 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127894 | MULTIPOLYGON (((2056301.039 5802553.702, 20562... |
| 388204 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127900 | MULTIPOLYGON (((2054727.599 5803954.788, 20547... |
| 445455 | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Gravel or Rock | 56 | 56 | 56 | 56 | 16 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb2000515912 | MULTIPOLYGON (((2055082.956 5803253.022, 20550... |
In [90]:
LCDB[LCDB.intersects(largest_area.geometry)].drop(columns="EditDate").explore("Name_2018", legend=True)
✔️ 358 ms (2024-12-02T09:45:43/2024-12-02T09:45:44)
Out[90]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [100]:
LCDB.loc[LCDB.intersects(largest_area.geometry), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True})
✔️ 300 ms (2024-12-02T09:51:28/2024-12-02T09:51:28)
Out[100]:
Class_2018 71 Wetland_18 False Onshore_18 True Name: 0, dtype: object
In [41]:
all_streams = REC.unary_union
✔️ 26.2 s (2024-12-02T08:49:13/2024-12-02T08:49:39)
In [57]:
features["distance_to_river"] = largest_area.geometry.distance(all_streams)
features["distance_to_river"]
✔️ 398 ms (2024-12-02T09:03:29/2024-12-02T09:03:29)
Out[57]:
0.0
In [113]:
def get_features(row):
geom = row.geometry
masked_old, transform = mask(old, [geom], nodata=np.nan)
nonan = np.argwhere(~np.isnan(masked_old[0]))
top_left = nonan.min(axis=0)
bottom_right = nonan.max(axis=0)
masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
masked_new, transform = mask(new, [geom], nodata=np.nan)
masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
raster = rasterize([geom], out_shape=result.shape, transform=new.transform)
masked_diff = np.where(raster, diff, np.nan)
masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
row = row.to_dict()
row.update({
"old_min": np.nanmin(masked_old),
"old_max": np.nanmax(masked_old),
"old_mean": np.nanmean(masked_old),
"old_median": np.nanmedian(masked_old),
"old_std": np.nanstd(masked_old),
"new_min": np.nanmin(masked_new),
"new_max": np.nanmax(masked_new),
"new_mean": np.nanmean(masked_new),
"new_median": np.nanmedian(masked_new),
"new_std": np.nanstd(masked_new),
"diff_min": np.nanmin(masked_diff),
"diff_max": np.nanmax(masked_diff),
"diff_mean": np.nanmean(masked_diff),
"diff_median": np.nanmedian(masked_diff),
"diff_std": np.nanstd(masked_diff),
"distance_to_river": geom.distance(all_streams)
})
row.update(LCDB.loc[LCDB.intersects(geom), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True}))
attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
attributes = xdem.terrain.get_terrain_attribute(
masked_diff,
resolution=old.res,
attribute=attribute_names
)
for i, name in enumerate(attribute_names):
row.update({
f"{name}_min": np.nanmin(attributes[i]),
f"{name}_max": np.nanmax(attributes[i]),
f"{name}_mean": np.nanmean(attributes[i]),
f"{name}_median": np.nanmedian(attributes[i]),
f"{name}_std": np.nanstd(attributes[i]),
})
return pd.Series(row)
#get_features(areas.iloc[0])
features = areas.progress_apply(get_features, axis=1)
features
✔️ 4 min 35 s (2024-12-02T10:03:06/2024-12-02T10:07:41)
0%| | 0/188 [00:00<?, ?it/s]
Out[113]:
| geometry | area | old_min | old_max | old_mean | old_median | old_std | new_min | new_max | new_mean | new_median | new_std | diff_min | diff_max | diff_mean | diff_median | diff_std | distance_to_river | Class_2018 | Wetland_18 | Onshore_18 | roughness_min | roughness_max | roughness_mean | roughness_median | roughness_std | slope_min | slope_max | slope_mean | slope_median | slope_std | aspect_min | aspect_max | aspect_mean | aspect_median | aspect_std | curvature_min | curvature_max | curvature_mean | curvature_median | curvature_std | terrain_ruggedness_index_min | terrain_ruggedness_index_max | terrain_ruggedness_index_mean | terrain_ruggedness_index_median | terrain_ruggedness_index_std | rugosity_min | rugosity_max | rugosity_mean | rugosity_median | rugosity_std | profile_curvature_min | profile_curvature_max | profile_curvature_mean | profile_curvature_median | profile_curvature_std | planform_curvature_min | planform_curvature_max | planform_curvature_mean | planform_curvature_median | planform_curvature_std | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 176 | POLYGON ((2054419.000 5803910.000, 2054419.000... | 530308.0 | 264.148010 | 520.244995 | 347.975189 | 350.695984 | 49.789028 | 264.705994 | 519.897034 | 356.345886 | 364.740967 | 53.155849 | -6.874023 | 43.521973 | 8.370939 | 4.911972 | 9.060170 | 0.000000 | 71 | False | True | 0.024994 | 16.997986 | 0.866442 | 0.501007 | 0.900825 | 0.000309 | 83.094493 | 15.825287 | 9.967977 | 14.601295 | 0.000000 | 359.999759 | 172.747728 | 183.009274 | 104.099068 | -1643.609619 | 2639.605713 | 1.731523 | 0.808716 | 43.725699 | 0.024669 | 27.121265 | 0.893986 | 0.540642 | 0.893833 | 1.000079 | 9.886722 | 1.108306 | 1.026156 | 0.199585 | -1767.249520 | 1433.146525 | -1.254866 | -0.584297 | 28.491703 | -624.432617 | 978.029107 | 0.476682 | 0.186515 | 21.372796 |
| 155 | POLYGON ((2054514.000 5804756.000, 2054514.000... | 319463.0 | 384.855011 | 949.697021 | 698.896851 | 715.611023 | 132.071274 | 384.348999 | 949.190002 | 684.987671 | 692.672974 | 129.693771 | -52.056030 | 5.609009 | -13.909487 | -9.703979 | 12.755947 | 0.000000 | 54 | False | True | 0.027039 | 10.799988 | 1.465731 | 1.246033 | 0.944709 | 0.021028 | 77.297658 | 26.350639 | 24.870094 | 14.383982 | 0.000000 | 359.999410 | 183.019811 | 185.773786 | 100.256831 | -954.620361 | 692.401123 | -1.137023 | -1.690674 | 58.625323 | 0.035931 | 10.752627 | 1.490565 | 1.282474 | 0.923966 | 1.000176 | 4.733187 | 1.209573 | 1.132283 | 0.229804 | -440.434936 | 585.141157 | 0.664059 | 0.707163 | 35.661075 | -369.479204 | 312.140260 | -0.473099 | -0.604222 | 31.241830 |
| 37 | POLYGON ((2056589.000 5808659.000, 2056589.000... | 304243.0 | 543.239014 | 968.459961 | 752.026978 | 758.137024 | 89.943939 | 542.731018 | 968.713989 | 750.921082 | 757.086975 | 90.210747 | -13.114990 | 5.141968 | -1.105809 | -0.944031 | 1.355662 | 0.000000 | 12 | False | True | 0.014038 | 6.578003 | 0.709351 | 0.576965 | 0.492890 | 0.019821 | 66.326472 | 13.694865 | 11.374013 | 9.495378 | 0.000000 | 359.998682 | 180.263642 | 180.474840 | 104.744639 | -472.998047 | 416.595459 | -1.971487 | -1.379395 | 39.778407 | 0.014390 | 5.968466 | 0.742228 | 0.609774 | 0.502908 | 1.000017 | 2.568007 | 1.063572 | 1.033145 | 0.087242 | -275.317190 | 316.405986 | 1.262352 | 0.660475 | 23.922942 | -335.581215 | 226.580603 | -0.709152 | -0.401965 | 21.456914 |
| 115 | POLYGON ((2054495.000 5806342.000, 2054495.000... | 173781.0 | 333.551025 | 461.656006 | 378.365387 | 375.186005 | 26.501810 | 334.056000 | 461.194000 | 381.083679 | 378.898010 | 26.786768 | -4.665009 | 12.177979 | 2.718257 | 2.453003 | 1.906919 | 0.000000 | 12 | False | True | 0.028015 | 6.747986 | 0.664326 | 0.454987 | 0.605888 | 0.028203 | 69.110675 | 12.514475 | 8.870377 | 11.030167 | 0.000000 | 359.999000 | 164.568461 | 177.295377 | 100.462405 | -553.802490 | 687.600708 | 2.607275 | 1.797485 | 41.001715 | 0.032552 | 7.412247 | 0.699980 | 0.489235 | 0.623502 | 1.000124 | 2.874868 | 1.065598 | 1.021682 | 0.119312 | -479.441262 | 366.685717 | -2.050307 | -1.421865 | 27.232389 | -357.169244 | 403.752029 | 0.556968 | 0.261994 | 19.498026 |
| 126 | POLYGON ((2054821.000 5806083.000, 2054821.000... | 147413.0 | 405.266998 | 984.378967 | 666.230042 | 680.250000 | 128.762848 | 404.644012 | 984.406982 | 662.956543 | 677.700012 | 128.024765 | -25.545959 | 10.151001 | -3.273446 | -2.276978 | 3.636282 | 0.000000 | 71 | False | True | 0.049988 | 8.898010 | 1.275577 | 1.083984 | 0.813165 | 0.040820 | 73.866107 | 23.143639 | 21.255515 | 13.418823 | 0.000000 | 359.998574 | 184.088882 | 176.004473 | 102.389142 | -767.303467 | 707.580566 | -4.001010 | -3.497314 | 62.809494 | 0.051240 | 8.480720 | 1.322321 | 1.132742 | 0.821058 | 1.000402 | 3.628650 | 1.171587 | 1.106380 | 0.191066 | -536.233840 | 418.001613 | 2.441112 | 1.776779 | 38.928714 | -393.931141 | 379.570854 | -1.559898 | -1.174689 | 33.025140 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 51 | POLYGON ((2057180.000 5807155.000, 2057180.000... | 4309.0 | 754.164001 | 862.838013 | 810.229065 | 812.798035 | 27.570267 | 753.013000 | 862.265991 | 808.988464 | 811.622009 | 27.611164 | -5.234009 | 0.836975 | -1.240615 | -1.079041 | 0.706411 | 377.665475 | 69 | False | True | 0.091980 | 3.499023 | 0.810893 | 0.724518 | 0.466198 | 0.258646 | 54.086741 | 15.535889 | 14.296695 | 9.054989 | 0.000000 | 359.873676 | 180.929910 | 181.396905 | 105.351676 | -295.800781 | 442.993164 | -7.090177 | -4.492188 | 47.193806 | 0.086022 | 3.618491 | 0.855655 | 0.758100 | 0.486536 | 1.001061 | 1.724303 | 1.077381 | 1.049964 | 0.088696 | -238.362703 | 175.115296 | 2.803598 | 1.571423 | 27.481605 | -176.705037 | 204.630461 | -4.286579 | -2.724701 | 26.676936 |
| 77 | POLYGON ((2057221.000 5806249.000, 2057221.000... | 4231.0 | 463.237000 | 558.276001 | 506.580353 | 501.701996 | 21.162735 | 462.401001 | 557.640991 | 504.420471 | 500.536987 | 21.546041 | -9.316010 | 0.420044 | -2.159984 | -1.524994 | 1.754438 | 0.000000 | 16 | False | True | 0.123962 | 5.431030 | 1.354514 | 1.208984 | 0.781550 | 0.243961 | 63.330647 | 24.899263 | 23.870917 | 12.951836 | 0.117864 | 359.932983 | 184.918589 | 175.719571 | 102.218871 | -365.005493 | 229.119873 | -11.565558 | -9.808350 | 54.859841 | 0.102363 | 4.997128 | 1.382453 | 1.254762 | 0.757155 | 1.000989 | 2.266732 | 1.179069 | 1.126430 | 0.172435 | -180.493919 | 211.376430 | 5.078110 | 3.214066 | 31.870435 | -207.633207 | 137.569661 | -6.487448 | -5.183216 | 31.373435 |
| 181 | POLYGON ((2054921.000 5802132.000, 2054921.000... | 4168.0 | 727.848999 | 777.617981 | 750.956482 | 751.011475 | 10.164408 | 727.346008 | 777.073975 | 748.487061 | 747.432495 | 10.804439 | -7.746033 | 0.098999 | -2.469438 | -1.833984 | 1.776451 | 313.839735 | 71 | False | True | 0.050049 | 2.488037 | 0.713715 | 0.641541 | 0.372780 | 0.161113 | 41.709637 | 14.652946 | 13.619203 | 7.531259 | 0.239895 | 359.927197 | 183.505383 | 172.413544 | 104.640565 | -209.301758 | 189.300537 | -3.036568 | -1.809692 | 25.295615 | 0.050731 | 2.201049 | 0.710312 | 0.651643 | 0.352425 | 1.000237 | 1.344907 | 1.051155 | 1.036017 | 0.048002 | -116.756886 | 137.333515 | 1.427443 | 0.626500 | 14.888978 | -90.982558 | 93.936800 | -1.609124 | -0.725124 | 13.837087 |
| 112 | POLYGON ((2056165.000 5805098.000, 2056165.000... | 4168.0 | 330.009003 | 414.976013 | 347.387390 | 338.135986 | 19.942204 | 329.310974 | 414.473022 | 345.676178 | 335.724487 | 20.337126 | -6.242035 | 0.404999 | -1.711227 | -1.103989 | 1.331259 | 20.061887 | 16 | False | True | 0.067993 | 5.648010 | 0.916036 | 0.627975 | 0.812676 | 0.072869 | 67.422201 | 16.667242 | 11.679819 | 14.067467 | 0.000371 | 359.989076 | 160.836444 | 161.757142 | 101.798624 | -645.401001 | 310.797119 | -10.793974 | -5.749512 | 54.868025 | 0.061296 | 7.057928 | 0.977414 | 0.664636 | 0.868992 | 1.000381 | 2.796020 | 1.118655 | 1.040430 | 0.184431 | -213.102706 | 411.064214 | 7.399746 | 3.237263 | 38.168409 | -249.806248 | 200.170430 | -3.394228 | -2.160090 | 26.729371 |
| 70 | POLYGON ((2056954.000 5806554.000, 2056954.000... | 4069.0 | 580.815002 | 632.517029 | 607.255920 | 607.020020 | 11.594285 | 580.307007 | 632.016968 | 605.696106 | 604.859009 | 11.600400 | -5.670044 | -0.108032 | -1.559806 | -1.331970 | 0.911266 | 44.852111 | 69 | False | True | 0.078003 | 2.404053 | 0.727969 | 0.664001 | 0.367533 | 0.283751 | 40.625595 | 14.404642 | 13.091004 | 7.748489 | 0.185128 | 359.981782 | 184.885313 | 190.514454 | 104.413637 | -207.397461 | 179.193115 | -6.322378 | -4.388428 | 34.206604 | 0.078687 | 2.184222 | 0.748892 | 0.687885 | 0.359458 | 1.000577 | 1.326722 | 1.056693 | 1.041069 | 0.049973 | -87.630775 | 112.695939 | 2.519276 | 1.309808 | 19.584720 | -115.848622 | 91.562340 | -3.803102 | -2.074324 | 19.314333 |
188 rows × 61 columns
In [114]:
features.crs = new.crs
✔️ 15.6 ms (2024-12-02T10:07:50/2024-12-02T10:07:50)
In [116]:
features.explore("Class_2018", legend=True)
✔️ 1.6 s (2024-12-02T10:07:59/2024-12-02T10:08:00)
Out[116]:
Make this Notebook Trusted to load map: File -> Trust Notebook